#include "VariantStats.h"
#include "Variant.h"
#include <cmath>
#include <iomanip>

using namespace std;
using namespace vcf;

bool isEqual(double dX, double dY) {
  const double dEpsilon = 0.000001; // or some other small number
  return fabs(dX - dY) <= dEpsilon * fabs(dX);
}

main(int argc, char *argv[]) {
  if (argc < 4) {
    cout << "invalid usage" << endl;
    return 1;
  }
  int bp_filter_start = atoi(argv[2]);
  int bp_filter_end = atoi(argv[3]);
  int pseudocount = 0;

  if (argc == 5) {
    pseudocount = atoi(argv[4]);
  }

  string chromosome;
  if (argc == 6) {
    chromosome = argv[5];
  }
  string vcf_file(argv[1]);
  map<int, int> p_distribution;
  VariantCallFile variantFile;
  variantFile.open(vcf_file);
  if (!variantFile.is_open()) {
    cerr << "Could not open file " << vcf_file << endl;
    return 1;
  }
  Variant var(variantFile);
  
  string chr;
  int bp;
  while(variantFile.getNextVariantNoParse(var, chr, bp)) {
    if (chr != chromosome) {
      continue;
    }
    var.parse(variantFile.line);
    if (var.position >= bp_filter_start && var.position <= bp_filter_end) {
      continue;
    }
    if (var.id[0] == 'r' && var.id[1] == 's') {
      map<string, map<string, vector<string> > >::iterator it;
      // Iterate over all samples, collecting information on the number
      // of hetero/homozygous samples, and the read counts of each
      // allele.
      int alt = 0, ref = 0;
      for (it = var.samples.begin(); it != var.samples.end(); ++it) {
	vector<string> genotypeInfo = it->second["GT"];
	assert(genotypeInfo.size() == 1);
	string genotype_string = genotypeInfo[0];
	if (genotype_string != "0/0") {
	  vector<string> readCounts = it->second["AD"];
	  assert(readCounts.size() == 2);
	  if (genotype_string[0] == '1') {
	    alt += atoi(readCounts[0].c_str());
	  } else {
	    ref += atoi(readCounts[0].c_str());
	  }
	  alt += atoi(readCounts[1].c_str());
	}
      }
      if (alt != 0) {
	int alt_promotion = (alt / ((float)(alt + ref))) * 100;
	p_distribution[alt_promotion] += 1;
      }
    }
  }

  
  map<int, int>::iterator iter;

  // Add pseudocounts
  for (iter = p_distribution.begin(); iter != p_distribution.end();
       ++iter) {
    pair<int, int> p = *iter;
    p_distribution[p.first] += pseudocount;
    cout << p.first << ":" << p_distribution[p.first] << endl;
  }

  int total_calls = 0;

  // Now total up calls
  for (iter = p_distribution.begin(); iter != p_distribution.end();
       ++iter) {
    pair<int, int> p = *iter;
    total_calls += p.second;
  }
  cout << "Total calls: " << total_calls << endl;
  // for (float f = 0.0; f <= 1; f += .01) {
  //   cout << f << ": " << p_distribution[f] / ((float)total_calls) << endl;
  // }
  int last = -1;
  for (iter = p_distribution.begin(); iter != p_distribution.end();
       ++iter) {
    pair<int, int> p = *iter;
    
    if (p.first - last != 1) {
      int j;
      for(j = last + 1; j < p.first; ++j) {
	cout << j << ": " << pseudocount / ((float)total_calls) << endl;
      }
    }
    cout << p.first << ": " << p.second/ ((float)total_calls) << endl;
    last = p.first;
  }
}
